"""
线性规划优化模型

问题描述：
    线性规划是一种在线性约束条件下优化线性目标函数的数学方法。
    标准形式：
    最大化(或最小化) c^T * x
    约束条件: A * x <= b, x >= 0
    其中c是目标函数系数向量，A是约束系数矩阵，b是约束右端向量

算法原理：
    线性规划使用单纯形法或内点法求解：
    1. 单纯形法：从可行域的一个顶点开始，沿着边移动到相邻顶点，
       直到找到最优解或确定问题无界
    2. 内点法：从可行域内部开始，通过势函数逼近最优解

核心概念：
    - 可行域：满足所有约束条件的解的集合
    - 最优解：使目标函数达到最优值的可行解
    - 对偶理论：每个线性规划问题都有对应的对偶问题
    - 敏感性分析：分析参数变化对最优解的影响

应用场景：
    线性规划广泛应用于：
    - 生产计划优化
    - 运输问题
    - 投资组合优化
    - 资源分配
    - 饮食搭配问题

作者：斯黄
日期：2025年3月4日
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linprog
from scipy.spatial import ConvexHull
import matplotlib.patches as patches

# 示例问题：生产计划优化
"""
某工厂生产两种产品A和B：
- 产品A每件利润：40元，需要原材料1单位，人工2小时
- 产品B每件利润：30元，需要原材料2单位，人工1小时
- 原材料供应：最多100单位
- 人工时间：最多80小时
- 求最大利润的生产计划
"""

print("=== 生产计划优化问题 ===")
print("目标：最大化利润")
print("产品A: 利润40元/件, 原材料1单位/件, 人工2小时/件")
print("产品B: 利润30元/件, 原材料2单位/件, 人工1小时/件") 
print("约束：原材料≤100单位, 人工≤80小时")
print()

# 定义线性规划问题
# 目标函数：max 40*x1 + 30*x2 (转换为min问题：min -40*x1 - 30*x2)
c = [-40, -30]  # 目标函数系数（取负数转换为最小化问题）

# 不等式约束 A_ub * x <= b_ub
# 1*x1 + 2*x2 <= 100  (原材料约束)
# 2*x1 + 1*x2 <= 80   (人工约束)
A_ub = [[1, 2],
        [2, 1]]
b_ub = [100, 80]

# 变量边界 (x1 >= 0, x2 >= 0)
x1_bounds = (0, None)
x2_bounds = (0, None)
bounds = [x1_bounds, x2_bounds]

# 求解线性规划
result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')

if result.success:
    print("=== 优化结果 ===")
    print(f"最优解: x1 = {result.x[0]:.2f}, x2 = {result.x[1]:.2f}")
    print(f"最大利润: {-result.fun:.2f} 元")
    print(f"原材料使用: {A_ub[0][0]*result.x[0] + A_ub[0][1]*result.x[1]:.2f}/100 单位")
    print(f"人工使用: {A_ub[1][0]*result.x[0] + A_ub[1][1]*result.x[1]:.2f}/80 小时")
else:
    print("优化失败:", result.message)

# 可视化可行域和最优解
def plot_linear_programming():
    plt.figure(figsize=(12, 8))
    
    # 定义坐标范围
    x1 = np.linspace(0, 60, 400)
    x2 = np.linspace(0, 60, 400)
    X1, X2 = np.meshgrid(x1, x2)
    
    # 约束条件
    # 原材料约束: x1 + 2*x2 <= 100
    constraint1 = (100 - X1) / 2
    # 人工约束: 2*x1 + x2 <= 80  
    constraint2 = 80 - 2*X1
    
    # 绘制约束线
    x1_line = np.linspace(0, 100, 100)
    plt.plot(x1_line, (100 - x1_line) / 2, 'r-', linewidth=2, label='原材料约束: x₁ + 2x₂ ≤ 100')
    plt.plot(x1_line, 80 - 2*x1_line, 'b-', linewidth=2, label='人工约束: 2x₁ + x₂ ≤ 80')
    
    # 可行域顶点
    vertices = np.array([[0, 0], [0, 50], [20, 40], [40, 0]])
    hull = ConvexHull(vertices)
    
    # 填充可行域
    for simplex in hull.simplices:
        plt.fill(vertices[simplex, 0], vertices[simplex, 1], alpha=0.3, color='lightblue')
    
    # 标记顶点
    for i, vertex in enumerate(vertices):
        plt.plot(vertex[0], vertex[1], 'ko', markersize=8)
        plt.annotate(f'({vertex[0]},{vertex[1]})', 
                    (vertex[0], vertex[1]), 
                    xytext=(5, 5), textcoords='offset points')
    
    # 绘制目标函数等值线
    profit_levels = [800, 1200, 1600, 2000]
    for profit in profit_levels:
        x2_profit = (profit - 40*x1_line) / 30
        plt.plot(x1_line, x2_profit, '--', alpha=0.7, label=f'利润 = {profit}')
    
    # 标记最优解
    if result.success:
        plt.plot(result.x[0], result.x[1], 'r*', markersize=15, label=f'最优解 ({result.x[0]:.1f}, {result.x[1]:.1f})')
    
    plt.xlim(0, 60)
    plt.ylim(0, 60)
    plt.xlabel('产品A数量 (x₁)')
    plt.ylabel('产品B数量 (x₂)')
    plt.title('线性规划问题可视化')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()

plot_linear_programming()

# 敏感性分析：分析原材料供应变化对最优解的影响
def sensitivity_analysis():
    print("\n=== 敏感性分析 ===")
    print("分析原材料供应量变化对最优利润的影响:")
    
    material_supplies = range(80, 121, 5)
    optimal_profits = []
    optimal_x1 = []
    optimal_x2 = []
    
    for supply in material_supplies:
        # 更新约束
        b_ub_new = [supply, 80]
        result_new = linprog(c, A_ub=A_ub, b_ub=b_ub_new, bounds=bounds, method='highs')
        
        if result_new.success:
            optimal_profits.append(-result_new.fun)
            optimal_x1.append(result_new.x[0])
            optimal_x2.append(result_new.x[1])
            print(f"原材料{supply}单位: 最优利润{-result_new.fun:.2f}元, "
                  f"x1={result_new.x[0]:.2f}, x2={result_new.x[1]:.2f}")
        else:
            optimal_profits.append(0)
            optimal_x1.append(0)
            optimal_x2.append(0)
    
    # 绘制敏感性分析结果
    plt.figure(figsize=(15, 5))
    
    plt.subplot(1, 3, 1)
    plt.plot(material_supplies, optimal_profits, 'b-o')
    plt.xlabel('原材料供应量')
    plt.ylabel('最优利润')
    plt.title('利润 vs 原材料供应量')
    plt.grid(True)
    
    plt.subplot(1, 3, 2)
    plt.plot(material_supplies, optimal_x1, 'r-o', label='产品A')
    plt.plot(material_supplies, optimal_x2, 'g-o', label='产品B')
    plt.xlabel('原材料供应量')
    plt.ylabel('最优生产量')
    plt.title('最优生产计划 vs 原材料供应量')
    plt.legend()
    plt.grid(True)
    
    plt.subplot(1, 3, 3)
    profit_increase = np.diff(optimal_profits)
    supply_increase = np.diff(material_supplies)
    marginal_profit = profit_increase / supply_increase
    plt.plot(material_supplies[1:], marginal_profit, 'purple', marker='s')
    plt.xlabel('原材料供应量')
    plt.ylabel('边际利润')
    plt.title('边际利润 vs 原材料供应量')
    plt.grid(True)
    
    plt.tight_layout()
    plt.show()

sensitivity_analysis()

# 对偶问题分析
def dual_problem_analysis():
    print("\n=== 对偶问题分析 ===")
    print("原问题：最大化生产利润")
    print("对偶问题：最小化资源价值")
    
    # 构造对偶问题
    # 对偶变量：y1(原材料影子价格), y2(人工影子价格)
    # 目标：min 100*y1 + 80*y2
    # 约束：1*y1 + 2*y2 >= 40, 2*y1 + 1*y2 >= 30, y1,y2 >= 0
    
    c_dual = [100, 80]  # 对偶目标函数系数
    A_ub_dual = [[-1, -2], [-2, -1]]  # 对偶约束矩阵(转换为<=形式)
    b_ub_dual = [-40, -30]  # 对偶约束右端(转换为<=形式)
    bounds_dual = [(0, None), (0, None)]
    
    result_dual = linprog(c_dual, A_ub=A_ub_dual, b_ub=b_ub_dual, bounds=bounds_dual, method='highs')
    
    if result_dual.success:
        print(f"对偶最优解: y1 = {result_dual.x[0]:.2f}, y2 = {result_dual.x[1]:.2f}")
        print(f"对偶最优值: {result_dual.fun:.2f}")
        print(f"原问题最优值: {-result.fun:.2f}")
        print(f"强对偶性验证: |原最优值 - 对偶最优值| = {abs(-result.fun - result_dual.fun):.6f}")
        print("\n影子价格解释:")
        print(f"原材料影子价格: {result_dual.x[0]:.2f} 元/单位")
        print(f"人工影子价格: {result_dual.x[1]:.2f} 元/小时")

dual_problem_analysis()

# 多目标线性规划示例
def multi_objective_example():
    print("\n=== 多目标优化示例 ===")
    print("目标1：最大化利润")
    print("目标2：最小化原材料使用")
    
    # 使用加权和法求解多目标问题
    weights = [(1, 0), (0.8, 0.2), (0.6, 0.4), (0.4, 0.6), (0.2, 0.8), (0, 1)]
    
    solutions = []
    for w1, w2 in weights:
        # 组合目标函数：w1*(-profit) + w2*(material_usage)
        # profit = 40*x1 + 30*x2, material = x1 + 2*x2
        c_multi = [-40*w1 + w2, -30*w1 + 2*w2]
        
        result_multi = linprog(c_multi, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')
        if result_multi.success:
            profit = 40*result_multi.x[0] + 30*result_multi.x[1]
            material = result_multi.x[0] + 2*result_multi.x[1]
            solutions.append((result_multi.x[0], result_multi.x[1], profit, material))
    
    # 绘制Pareto前沿
    profits = [sol[2] for sol in solutions]
    materials = [sol[3] for sol in solutions]
    
    plt.figure(figsize=(10, 6))
    plt.plot(materials, profits, 'ro-', linewidth=2, markersize=8)
    for i, (x1, x2, profit, material) in enumerate(solutions):
        plt.annotate(f'({x1:.1f},{x2:.1f})', (material, profit), 
                    xytext=(5, 5), textcoords='offset points')
    
    plt.xlabel('原材料使用量')
    plt.ylabel('利润')
    plt.title('多目标优化：Pareto前沿')
    plt.grid(True)
    plt.show()
    
    print("Pareto最优解集:")
    for i, (x1, x2, profit, material) in enumerate(solutions):
        print(f"解{i+1}: x1={x1:.2f}, x2={x2:.2f}, 利润={profit:.2f}, 原材料={material:.2f}")

multi_objective_example() 